//-------------------------------------------------------------------------//
//                                                                         //
//  This benchmark is an OpenMP C version of the NPB SP code. This OpenMP  //
//  C version is developed by the Center for Manycore Programming at Seoul //
//  National University and derived from the OpenMP Fortran versions in    //
//  "NPB3.3-OMP" developed by NAS.                                         //
//                                                                         //
//  Permission to use, copy, distribute and modify this software for any   //
//  purpose with or without fee is hereby granted. This software is        //
//  provided "as is" without express or implied warranty.                  //
//                                                                         //
//  Information on NPB 3.3, including the technical report, the original   //
//  specifications, source code, results and information on how to submit  //
//  new results, is available at:                                          //
//                                                                         //
//           http://www.nas.nasa.gov/Software/NPB/                         //
//                                                                         //
//  Send comments or suggestions for this OpenMP C version to              //
//  cmp@aces.snu.ac.kr                                                     //
//                                                                         //
//          Center for Manycore Programming                                //
//          School of Computer Science and Engineering                     //
//          Seoul National University                                      //
//          Seoul 151-744, Korea                                           //
//                                                                         //
//          E-mail:  cmp@aces.snu.ac.kr                                    //
//                                                                         //
//-------------------------------------------------------------------------//

//-------------------------------------------------------------------------//
// Authors: Sangmin Seo, Jungwon Kim, Jun Lee, Jeongho Nah, Gangwon Jo,    //
//          and Jaejin Lee                                                 //
//-------------------------------------------------------------------------//

#include "header.h"
#include "mpinpb.h"
#include <math.h>

//---------------------------------------------------------------------
// this function computes the norm of the difference between the
// computed solution and the exact solution
//---------------------------------------------------------------------
void error_norm(double rms[5])
{
  int i, j, k, m, d;
  double xi, eta, zeta, u_exact[5], add;
  double rms_local[5];
  int scount, rcount,kstart, kend,kgap, krem;
  for (m = 0; m < 5; m++) {
    rms[m] = 0.0;
  }

 // #pragma omp parallel default(shared) \
   //       private(i,j,k,m,zeta,eta,xi,add,u_exact,rms_local) shared(rms)
 // {
    for (m = 0; m < 5; m++) {
      rms_local[m] = 0.0;
    }

   // #pragma omp for nowait
  kgap= grid_points[2]/ total_nodes;
  krem = grid_points[2]% total_nodes;
  if(krem == 0){
	kstart=node*kgap;
	kend= kstart+kgap;}
  else {
	  if(node+1 <= krem) 
	  {
		kstart = node*(kgap+1); //1 , 4
		kend= kstart+kgap+1; // 4, 7,
	  }
	  else
	  {
		kstart= krem*(kgap+1)+(node - krem)*kgap; //	7 , 9
		kend = kstart +kgap; // 9, 11
	  }
	}
  for (k = kstart; k < kend; k++) {
     zeta = (double)k * dnzm1;
       for (j = 0; j <= grid_points[1]-1; j++) {
        eta = (double)j * dnym1;
	 for (i = 0; i <= grid_points[0]-1; i++) {
          xi = (double)i * dnxm1;
          exact_solution(xi, eta, zeta, u_exact);
          for (m = 0; m < 5; m++) {
            add = u[k][j][i][m]-u_exact[m];
            rms_local[m] = rms_local[m] + add*add;
	   //printf("add : %0.20f rms_local %0.20f\t",add, rms_local[m]);	
 	    
          }
        }
      }
    }
 
  //MPI_Barrier(MPI_COMM_WORLD);
  MPI_Allreduce(rms_local,rms, 5, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  
   /*for (m = 0; m < 5; m++) {
   printf("rms %0.20f\n", rms[m]);
     // #pragma omp atomic
      //rms[m] += rms_local[m];
    }*/
 // } //end parallel

  for (m = 0; m < 5; m++) {
    for (d = 0; d < 3; d++) {
      rms[m] = rms[m] / (double)(grid_points[d]-2);
    }
    rms[m] = sqrt(rms[m]);
  }
}


void rhs_norm(double rms[5])
{
  int i, j, k, d, m;
  double add;
  double rms_local[5];
   int scount, rcount,kstart, kend,kgap, krem;
  for (m = 0; m < 5; m++) {
    rms[m] = 0.0;
  }

  //#pragma omp parallel default(shared) private(i,j,k,m,add,rms_local) \
                                       shared(rms)
  //{
    for (m = 0; m < 5; m++) {
      rms_local[m] = 0.0;
    }
  kgap= grid_points[2]/ total_nodes;
  krem = grid_points[2]% total_nodes;
  if(krem == 0){
	kstart=node*kgap;
	kend= kstart+kgap;}
  else {
	  if(node+1 <= krem) 
	  {
		kstart = node*(kgap+1); //1 , 4
		kend= kstart+kgap+1; // 4, 7,
	  }
	  else
	  {
		kstart= krem*(kgap+1)+(node - krem)*kgap; //	7 , 9
		kend = kstart +kgap; // 9, 11
	  }
	}
   if(node == root)
	kstart = kstart+1;
   if(node == total_nodes-1)
	kend = kend-1;	
    //#pragma omp for nowait
    for (k = kstart; k < kend; k++) { 
      for (j = 1; j <= ny2; j++) {
        for (i = 1; i <= nx2; i++) {
          for (m = 0; m < 5; m++) {
            add = rhs[k][j][i][m];
            rms_local[m] = rms_local[m] + add*add;
          } 
        } 
      } 
    } 
  MPI_Allreduce(rms_local,rms, 5, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
    /*for (m = 0; m < 5; m++) {
      #pragma omp atomic
      rms[m] += rms_local[m];
    }
  } //end parallel*/

  for (m = 0; m < 5; m++) {
    for (d = 0; d < 3; d++) {
      rms[m] = rms[m] / (double)(grid_points[d]-2);
    }
    rms[m] = sqrt(rms[m]);
  }
}

